Comparison of Dissipative Particle Dynamics and Langevin thermostats 
for out-of-equilibrium simulations of polymeric systems 
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In this work we compare and characterize the behavior of Langevin and Dissipative Particle Dy- 
namics (DPD) thermostats in a broad range of non-equilibrium simulations of polymeric systems. 
Polymer brushes in relative sliding motion, polymeric liquids in Poiseuille and Couette flows, and 
brush-melt interfaces are used as model systems to analyze the efficiency and limitations of different 
Langevin and DPD thermostat implementations. Widely used coarse-grained bead-spring models 
under good and poor solvent conditions are employed to assess the effects of the thermostats. We 
considered equilibrium, transient, and steady state examples for testing the ability of the ther- 
mostats to maintain constant temperature and to reproduce the underlying physical phenomena in 
non-equilibrium situations. The common practice of switching-off the Langevin thermostat in the 
flow direction is also critically revisited. The efficiency of different weight functions for the DPD 
thermostat is quantitatively analyzed as a function of the solvent quality and the non-equilibrium 
situation. 



I. INTRODUCTION 

The need and use of thermostats in computer simula- 
tions started with the beginning of the field itself. The 
original Molecular Dynamics (MD) method, intended for 
microcanonical ensemble simulations, was soon extended 
to different ensembles in order to mimic conditions in 
which experiments are actually performed. The thermo- 
stat in MD simulations implies the assumption that the 
system transports heat "instantaneously" fast on the spa- 
tial scale of the simulation. Even when this is arguably 
not completely correct in a real system many studies 
are faced with the situation of performing simulations at 
constant temperature as a way of obtaining a physically 
meaningful condition. 

It is a non-trivial challenge to achieve a constant tem- 
perature in a simulation of driven soft matter systems, 
such as polymer brushes interacting with flowing poly- 
mer melts, which we consider in the present work. 

In this work, we present simulation data which we ob- 
tained by integrating Langcvin-likc equations with the 
standard MD integrators using cither a Langevin or a 
DPD thermostat to maintain temperature. This is often 
called Stochastic Dynamics (SD) or Brownian Dynamics 
(BD)ii2i2, but we address a physical regime in which the 
friction and stochastic forces, added to the conservative 
forces of the system, are small enough in order not to sig- 
nificantly perturb the natural dynamics of the polymeric 
system. This is typically achieved by using the smallest 
possible value for the friction constant, 7, providing that 
the temperature is conserved under the desired physical 
conditions. This approach also requires that the system 
has intermediate to high monomer densities to warrant 
that the friction due to the conservative interactions is 
significantly larger than the frictional and random forces 
arising from the thermostat. The regime of dilute poly- 



mer solutions is excluded from this study because, in this 
case, the polymer-solvent interactions are important for 
the local dynamics of the molecules. 

In this sense we consider the Langevin and DPD fric- 
tional and stochastic forces as thermostats that are added 
to the true conservative forces of the system, as it is gen- 
erally employed in MD. This physical regime is of great 
interest for a number of soft matter systems, such as poly- 
meric interfaces, blends and melts, that are successfully 
studied in the framework of coarse-grained models^££ 

A case which has attracted abiding interest is the sim- 
ulation of out-of-equilibrium phenomena in which a rate 
of energy must be injected into the system to drive it out 
of equilibrium. This energy must be removed in order to 
keep the temperature constant, and this is usually done 
by the action of a thermostat. There are excellent reviews 
that describe the different types of thermostats and their 
respective advantages and limitations for studying vari- 
ous systems and physical phenomenapii^ Of course, it is 
crucial that the dynamical behavior observed in a sim- 
ulation faithfully represents the actual dynamics of the 
desired system and is essentially free from artifacts in- 
troduced by the thermostat. This issue shall be explored 
for out-of-equilibrium simulations of polymeric systems 
in our study. 

In this article, we focus on the behavior of Langevin 
and DPD^ i 10 ! 11 ! 12 thermostats for a range of typical poly- 
meric systems in non-equilibrium conditions. The former 
has been widely used in equilibrium simulations but is 
known to have undesirable properties, such as screen- 
ing of hydrodynamic interactions and lack of Galilean 
invariancejii^ in non-equilibrium situations. A typical 
workaround when using the Langevin thermostat for non- 
equilibrium simulations consists in switching-off the ther- 
mostat in the direction in which non-conservative exter- 
nal forces are applied to the system or applying it only in 
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one spatial direction. In this way one recovers momen- 
tum conservation in the shear direction, while conserving 
the temperature by applying the thermostat in the per- 
pendicular direction in which no direct non-conservative 
force is applied££^ 

The DPD scheme only recently has started to be uti- 
lized as a standalone thermostatic It was originally de- 
veloped as a method for performing meso-scale simula- 
tions by combining this thermostat and very "soft" po- 
tentials. The latter allow for the use of a large time step 
in MD simulations.— i 10 ' 11 The maximal time step that 
is permissible in DPD simulations has been investigated 
thoroughlyi 8 i 16 i 17 i 18 Utilizing a DPD thermostat in con- 
junction with "hard" potentials - typical of coarse-grained 
models widely used for polymers and other condensed 
matter systems 1 ^ - one looses this advantage, and one 
must take a time step on the order of that typically used 
in MD simulations of coarse-grained models. The local 
conservation of momentum and the Galilean invariance, 
however, are inherited from the original DPD method, 
and possibly this is a great advantage. 

In the following, we consider three polymeric systems 
to assess the effects of two versions of Langevin ther- 
mostats (with and without switching-off the thermostat 
in the flow direction) and the DPD thermostat: (a) single 
end-grafted polymer layers (so-called polymer brushes), 
(b) two opposing and intcrdigitating polymer brushes, 
and (c) a brush-melt interface, which exhibits a rich wet- 
ting behavior. The equilibrium properties of these ref- 
erence systems are interesting and have been compre- 
hensively studied in previous work a 4 ' 15 i 20 ' 21 i 22 i 23 ' 24 . The 
brush layers are characterized by the number of grafted 
chains per unit area or grafting density p g . Typical 
equilibrium density profiles for these three systems are 
shown in Fig. [TJ Additionally, a simple bulk system of 
a polymeric liquid with periodic boundary conditions in 
all three directions was considered. 

These systems are also of great interest out of equi- 
librium, for instance, because of the surprisingly small 
friction of two opposing brushes sliding past each 
others.- 21 - 2 ^ 

The interface between a brush and a melt of identical 
chains is a prototypical example for a copolymer-laden 
interface or a melt in contact with a soft, elastically dc- 
formablc substrate (e.g., confining brush-coated walls of 
a channel).— In addition to the rich wetting properties, 
typical applications (e.g., droplet break-up in a polymer 
blend under shear or flow in a microfluidic channel) in- 
volve flow and shear at the brush-melt interface. There- 
fore, the study of boundary conditions and the rheological 
properties of the macromolecular liquid subjected to dif- 
ferent types of flows make the non-equilibrium properties 
of this system particularly interesting. 

The behavior of the thermostats in equilibrium, tran- 
sient, and different kinds of steady states was tested in 
a wide number of typical situations that can be encoun- 
tered in simulations. Also, Poiseuille and Couette flows 
of the polymeric liquid were considered to compare dif- 



ferent weight functions of the DPD thermostat. 

The details of the thermostats and the polymer model 
are explained in section HH Section IIIII presents the dis- 
cussion of our results, which begins with a quantitative 
study of the relative strengths of thermostats for a given 
set of parameters. This section is divided in subsections 
corresponding to each different system: the analysis of 
single-brush layer transient states is presented in section 
IIII A[ and the steady state of two polymer brushes in 
relative sliding motion is discussed in IIII Bl for two dif- 
ferent grafting densities. In this way, we address two 
regimes: concentrated solution or melt in which hydro- 
dynamic interactions are screened and the a more dilute 
regime. Finally, the study of different DPD weight func- 
tions and their efficiency for conserving the temperature 
in the strong out-of-cquilibrium regime for Couette and 
Poiseuille flows of the brush-melt interface is described 
in subsection IIII CI Discussion and concluding remarks 
are presented in section ITVl 



II. 



POLYMER MODEL AND THERMOSTAT 
DETAILS 



We used a well established coarse-grained bead-spring 
models^ for polymers with excluded volume and intra- 
molecular interactions. This model has been applied to a 
variety of thermodynamic conditions, chain lengths, and 
physical regimes such as glasses, melts, dilute solutions, 
etci^ ' 27 i 28 ' 29 The interaction between neighboring beads 
along the same polymer is modeled by a finite extensible 
non-linear elastic (FENE) potential: 
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where the maximum allowed bond length is Rq = 1.5(7, 
the spring constant is k = 30e/<7 2 , and r = |r.; — r j | de- 
notes the distance between neighboring monomers. Ex- 
cluded volume interactions at short distances and van- 
der-Waals attractions between segments are described by 
a truncated and shifted Lennard-Jones (LJ) potential: 



with 



U(r) = U hJ (r)-Uu(r c ), 



U L] {r)=As 



(2) 



(3) 



where the LJ parameters, £ = 1 and a = 1, define the 
units of energy and length, respectively. The temper- 
ature is therefore given in units of e/ks, with ks the 
Boltzmann constant. U^j(r c ) is the LJ potential eval- 
uated at the cut-off radius. We considered two values 
as cut-off distance: (i) twice the minimum of the LJ 
potential: r c = 2 x 2« ~ 2.24er, which allows to con- 
sider poor solvent conditions, and (ii) r c = 2s ~ 1.12cr, 
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which models good solvent conditions. In the latter case, 
the interactions between monomers of different chains are 
purely repulsive , 21 i 22 i 23 whereas in the former case, longer 
ranged attractions are included giving rise to liquid- vapor 
phase separation and droplet formatio n 30 ' 31 below the 0- 
temperature, = 3.3e/fcB- We analyze the efficiency of 
the thermostats for both cases. 

The substrate is modeled as an idealized flat and im- 
penetrable wall, which interacts with the polymer seg- 
ments via an integrated Lennard-Jones potential : 

W*) = 14(f) 9 - A (£) 3 ' ( 4 ) 

where A = 3.2e, used throughout this work, is sufficient 
to make the liquid wet the bare substrate. 31 ' 32 The teth- 
ered beads are fixed randomly in the grafting plane at a 
distance of 1.2c from the wall position for all the cases. 
In the following, we will use LJ units^ for all quantities, 
unless explicitly mentioned otherwise. 

A. Langevin and Dissipative Particle Dynamics 
thermostats 

Both, Langevin and DPD thermostats, can be written 
in a general form, starting from the Hamiltonian equa- 
tions of motion^ 

Pi 

= — 

rrii 

Pi = F. + RD + F*, (5) 

where Fj is the total conservative force on each parti- 
cle and rrii the mass of each particle. Fj D and Fi R are 
the forces due to the thermostat and will be of differ- 
ent form for Langevin and DPD cases. The difference 
between DPD and Langevin thermostats is the way in 
which random and dissipative forces are applied. 

In the case of Langevin thermostats, the dissipative 
force on particle i is given by F; D = — 7V,, where 7 is 
the friction coefficient and the particle velocity. The 
random force, Ff-, has zero mean value and its variance 
satisfies 3 - 

(F*(t)F*(t')) = 2 1 Tk B 8 ij 8 liV 8(t - f) , (6) 

where the indices i and j label particles, fj, and v Carte- 
sian components, and T is the temperature at which the 
system is simulated. 

For the DPD case^tL - the dissipative and frictional 
forces are applied in a pair-wise form, such that the sum 
of thermostating forces acting on a particle pair equals 
zero. The expression for the forces is the following: 

= X I'';.; ■ I"'!! -7" D (rtf)(*« ■ v )f - 
Fi R = £F* ; Fg=^(ry)<?«% (7) 



where for each vector a we define ay = a; — a^, 7 is 
the friction constant, and a the noise strength. Friction 
and noise, 7 and a, obey the relation a 2 = 2ksTj, and 
the associated weight functions satisfy the fluctuation- 
dissipation theorem if the following relation is fulfilled: 11 

= (8) 

8ij is a random variable with zero mean and second mo- 
ment 

(Oij(t)e kl (t)) = (SijSji + 5aS jk )S{t - t). (9) 
The standard weight functions found in the literature are: 

[^f = ^ = ^- r l r ^ r<r < < (10) 

[0, r > r c 

where r c is the cut-off radius for a given molecular model. 
However, we emphasize that Eq. (jTU)) is just the typical 
choice when the DPD thermostat is employed in conjunc- 
tion with "soft" potentials. For arbitrary models, one can 
choose a different set of functions providing that they ful- 
fill Eq. ([5]) , and one can exploit this freedom to optimize 
the efficiency of the thermostat for "hard" potentials. In 
this work, we will use the standard weights, but also test 
other possibilities, whose forms are given in the first row 
of table 12 The equations of motion [Eq. were inte- 
grated using the velocity Verlet algorithmic with a time 
step of dt = 0.002t, where r = a(mfe) 1 ^ 2 denotes the 
time unit in terms of LJ parameters. 

III. RESULTS 

In this section, the results corresponding to the three 
polymeric systems, whose equilibrium density profiles are 
shown in Fig.[TJ will be analyzed. The difference between 
weight functions and the way random and friction forces 
are applied in Langevin and DPD thermostats does not 
allow for a direct comparison of the friction strength 7 
between both schemes. To obtain a direct measurement 
of thermostat strengths, we computed the mean friction 
and dissipative forces as a function of 7 for a polymeric 
liquid of 10-bead chains in a bulk solution (using periodic 
boundary conditions in all spatial directions) . In the case 
of the DPD thermostat, the standard weight functions 
were used [see Eq. (fTD])], Figure [2] shows the total force as 
a function of 7. Poor and good solvent conditions exhibit 
quite a different behavior. In the first case, the Langevin 
and DPD thermostat show a very similar behavior for the 
whole range of 7. For good solvent conditions (only the 
repulsive part of the LJ potential is kept), however, the 
mean Langevin forces are two orders of magnitude larger 
than those of the corresponding DPD counterpart (for 
the standard choice of weight functions). This means, for 
example, that a friction constant 7dpd = 2 is equivalent 
to a value of 7lgv = 0.01, as regards the mean value of 
the "thermostat force" acting on each bead. The ratio 
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Figure 1: (Color online) Equilibrium density profiles for the 
three systems under consideration: (a) single brush layer, (b) 
two opposing brush layers in interaction at two different graft- 
ing densities, and (c) brush-melt interfaces. The insets show 
the product of upper and lower brush density profiles p u x p\ 
(b) or the product of brush and melt phases pb x p m (c), ac- 
counting for the level of interdigitation of the different phases. 
In all cases, the temperature is T = 1.68. The chain length is 
N = 30 for cases (a) and (b) and N = 10 for case (c). The 
distance between the grafted beads of the opposing brushes is 
D = 17.5 for case (b) and D = 30 for case (c). The normal- 
ization is given by po — p g /D. 



(Fdpd) as a f unc ti° n °f 7 is shown with open symbols for 
sake of comparison. The reason for the big difference, 
in the case of good solvent, is the structure of the liquid 
and, more important, the small cut-off radius - there are 
very few beads in the range of the weight functions, and 
the standard weight functions are small in the vicinity 
of r c where most of the neighbors are located. The pair 
correlation functions and corresponding cut-off radii are 
shown in the inset of Fig. [2] 

At this point, it is important to recall that the origi- 
nal reason for choosing those weight functions [Eq. (fTO)) ] 
was based on the idea of using DPD together with "soft" 
potential s 9 ' 12 to achieve the largest possible time step. To 
this end, the thermostat forces need also to be smoothly 
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Figure 2: (Color online) Comparison of mean random and 
dissipative forces for DPD and Langevin (LGV) thermostats 
as a function of the friction constant, 7, for a bulk system with 
temperature T = 1.68 and density p = 0.61. The ratio |^gv) 
is shown in open circles for both, good solvent (gs) and poor 
solvent (ps) conditions. Inset: Pair correlation function <?2(r) 
of the polymeric liquid for good (solid line) and poor (dashed 
line) solvent conditions. The arrows indicate the position of 
the cut-off radius in each case. 



varying functions of the position in order to have the 
same properties as the conservative forces. Actually, the 
only constraint the weight functions must fulfill is the 
fluctuation-dissipation theorem**^ i.e., Eq. ([8]). 

We will see below that the standard choice can be even 
bad for non-equilibrium simulations, in which a signifi- 
cant amount of heat per unit time has to be removed. 

To quantify the efficiency of the thermostat to main- 
tain constant temperature, we define the number of ther- 
mostated particles, Atp, for a given DPD weight function 
pair as 

N TP =pof u R (r)g 2 (r)4irr 2 dr, (11) 
Jo 

where r c is the cut-off radius of the conservative poten- 
tials which coincides with the cut-off of the weight func- 
tions for the dissipative and frictional forces. g2 (r) is the 
pair correlation function for particles in the polymeric liq- 
uid. For the standard weight functions and good solvent 
conditions, Ntp is rather small (cf. Fig. [5]). 

The pair correlation functions for poor and good sol- 
vent conditions were taken from bulk simulations of a 
10-bead polymeric liquid at T = 1.68 and p = 0.61 which 
corresponds to the density of a melt that coexists with 
its vapor i^I In table U the calculated values of Atp, as 
given by Eq. (fTTj) . are shown for three different choices 
of weight functions. The second column shows Atp for 
the standard weight functions [Eq. (fTO)) ] . The third col- 
umn shows the square root of the usual weight functions 
which slightly increases the force in the region r < r c . 
The last column corresponds to constant weight func- 
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N tp (lj r ) 


u, R = 1 - X. 


V 


w R = 9(r c - r) 


good solvent 
(r c = 1.12) 


0.301 


0.904 


2.997 


poor solvent 
{r c = 2.24) 


6.753 


12.814 


28.868 



Table I: Number of thermostated particles, Ntp, for different 
weight functions using the DPD thermostat. See text and 
Eq. ([11]). 



tions, i.e., u> R = up = 0(r c — r), with 9 being Heav- 
iside's step function. The differences among the weight 
functions is evident from the values of Ntp- Different 
weight functions give rise to significant changes in the 
efficiency of the thermostat. In particular, the standard 
weight functions inherited from DPD models with "soft" 
potentials present the lowest value for A^tp which is sig- 
nificantly smaller than the values for other choices. As 
will be shown in section llll CI this is indeed an important 
issue in out-of-cquilibrium simulations. 



A. Single brush 

We consider a single polymer brush layer with TV = 30 
beads per chain in equilibrium under good solvent condi- 
tions. The brush stretches freely according to the balance 
of entropy and steric repulsion between the monomers A 
At time t = a constant wall velocity of v x = 1 is 
switched on in one direction, and the transient behavior 
of the brush is monitored. We took mean values over 10 
simulations starting from independent equilibrium con- 
figurations. For each simulation, 3 • 10 5 to 6 • 10 5 steps 
with a time step dt = 0.002 were performed. 

The rheological response of the polymer brush is ana- 
lyzed for three different cases: the usual Langevin ther- 
mostat with the same value of 7 in all spatial direc- 
tions, the Langevin thermostat with zero friction con- 
stant in shear direction (denoted as j x = 0), and the 
DPD thermostat. It is known that the usual Langevin 
thermostat (7^ ^ 0) does not reproduce the hydrody- 
namic behavior correctly because it does not conserve 
momentumiii^ and biases the flow profile in shear direc- 
tion. A common workaround to partially overcome these 
problems in non-equilibrium simulations of simple (lam- 
inar) flows is to switch-off the Langevin thermostat in 
the direction in which external non-conservative forces 
are applied.—^ This corresponds to our second approach 
with j x = 0. Alternatively, Langevin thermostat is fre- 
quently switched-off in two Cartesian components, being 
active only in the vorticity direction^^^ 

The thermostat's action can be rather understood as 
an implicit (ficticious) solvent acting on the polymer 
beads. In our coarse-grained simulation approach an ef- 
fect of the implicit solvent on the dynamical properties 
is undesirable. As we will show in the following, how- 
ever, there arc cases where such effects cannot be avoided. 



In the case of the Langevin thermostat, this solvent is 
at rest in the laboratory frame. When the brush layer 
starts moving through the implicit solvent, the behav- 
ior of Langevin and DPD thermostats differs drastically. 
A first evidence of these differences is shown in Fig. [3l 
where Langevin and DPD thermostats are compared us- 
ing a friction constant of 7 = 2 in all spatial directions. 
The angle a between the vector normal to the substrate 
and the average end-to-end vector, R c = ri — r_/v, with 
17 and r^r denoting the position vectors of the grafted 
and free end of a polymer chain, respectively, is shown 
as a function of time. While for the DPD thermostat 
a exhibits a decaying oscillatory behavior which ends in 
a steady state with the brush perpendicular to the wall 
(a = 0), the Langevin thermostat shows an angle which 
monotonously increases to a steady state value of 60°. 
This can be understood in terms of the lack of Galilean 
invariance of the Langevin thermostat. The brush is 
dragged through a ficticious solvent, which is always at 
rest in the laboratory frame. In ease of the DPD thermo- 
stat, the Galilean invariance, which follows from momen- 
tum conservation, implies that the polymer brush at rest 
and in steady state are equivalent. This case provides an 
example where the particular thermostat implementation 
plays a crucial role for the dynamics of the system and 
can even alter the results qualitatively. 

The oscillation frequency of a is a general property 
of the brush layer from which we can extract a char- 
acteristic response time. The component of the end- 
to-end vector in shear direction, R® (dashed lines in 
Fig. [3]) , and the mean shear stress (force per surface area) 
on the grafted heads of the polymer brush (see Fig. |4|) 
present a similar behavior. Defining for each observable 
A = A exp(— t/r c ), the decaying envelope of the oscillat- 
ing curves, maxima and minima can be brought onto the 
same curve (see Fig. [5]), yielding a characteristic time, 
t c = 53.94r. The j x = case for the Langevin thermo- 
stat (not shown) presents the same behavior as the DPD 
thermostat with the same characteristic time within the 
error. 

As observed for the mean thermostat forces (see Fig. [2]) 
in the good solvent case, similar forces for Langevin and 
DPD thermostats are obtained for values of 7 which differ 
by two orders of magnitude. We therefore performed the 
simulations for the Langevin thermostat with 7 = 0.01 
in all spacial directions, which corresponds to 7 = 2 for 
the DPD case. Indeed, we observe the oscillations also 
for the Langevin thermostat for a sufficiently small value 
of 7 (see Fig. HJ. This shows that the over-damping in 
the previous case (7 = 2) is due to the large, but not 
uncommon value of 7. 

While the oscillations are reproduced for the smaller 
value of 7, the steady state shear stress remains finite. It 
can easily be calculated for both values of 7 via 



a° s = Ps jNmv x , (12) 
with m = 1 the monomer mass. As expected, the steady 
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Figure 3: (Color online) Transient evolution of the inclination 
angle a of the polymer brush (solid line) and the x component 
of the end-to-end vector (dashed line) for Langevin and 
DPD thermostats. r c is the characteristic response time (see 
text and Fig. [SJ) 



state for j x — yields <r s — a® = 0, shown as a dashed 
line in the inset of Fig. 2J 

Figure [5] shows the transient behavior of the normal 
pressure exerted by the brush layer. Again, oscillations 
of the stress are observed for the DPD thermostat, while 
they are not directly observed in either version of the 
Langevin thermostat even for the smaller value of the 
friction coefficient. A Fourier frequency analysis of the 
time sequence, however, exhibits a peak at the same fre- 
quency for all cases. 

At first one could argue that the Langevin thermo- 
stat over-damps the oscillations because local momen- 
tum is only conserved in the shear direction, and a non- 
equilibrium situation in which there is a strong coupling 
between directions cannot be faithfully reproduced. This 
coupling would be reinforced by the chain connectivity of 
the polymeric system. We found, however, that for good 
solvent conditions and 7 = 0.5 the DPD thermostat is not 
able to maintain constant temperature with the standard 
choice of weight functions. Taking, for example, weight 



functions as 



,R _ 



= lu d = Q(r 



") (see fourth column of 



table H]), the DPD thermostat maintains the temperature 
for a wall velocity interval v x <G [0,1], but also the normal 
pressure oscillations are suppressed. 



B. Two opposing brushes 

A system of two opposing brush layers was studied un- 
der constant shear. It is similar to that already studied 
in previous workS ) 21 i 22 ' 23 i 34 ' 35 i 36 and originally attracted 
much interest because of its extraordinary small lateral 
friction forces. The equilibrium density profiles are de- 
picted in Fig. QJb) . We considered brush layers of chains 
with N — 30 beads at two different grafting densities: 
p g = 1.2/9* and p g = 4.9p*, where p* = l/vrfl 2 . (R g = 3.02 
being the radius of gyration of a single chain in so- 
lution) is the grafting density characterizing the grad- 
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Figure 4: (Color online) Comparison of shear stress for DPD 
and Langevin thermostats. For the second case two different 
values of 7 were considered. The inset shows the systematic 
difference for DPD and Langevin thermostat with 7 = 0.01, 
close to steady state, and the approach to <r s = for the 
Langevin thermostat with 7^ = (dashed line). r c is similar 
to Fig. [3] The mean stress for Langevin is indicated with a 
horizontal dashed lines to improve clarity. 
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Figure 5: (Color online) Maxima and minima of the oscil- 
lating curves (Figs. [3] and 2} for the single brush system. A 
characteristic frequency of the system is found when using 
DPD or Langevin thermostats with 7 = 0.01. 
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Figure 6: (Color online) Upper panel: Normal stress as a 
function of time for DPD. Lower panel: Langevin thermostat 
with 7 = 0.01 in all directions and -y x — along the shear 
direction. The curve corresponding to the latter case was 
shifted by 0.2 to improve clarity. 
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ual crossover from the mushroom to the brush regime. 
Pg = 1.2p* is a system within the crossover regime be- 
tween mushroom and brush, whereas the latter choice 
of p g leads to a semi-dilute brush^ For both consid- 
ered grafting densities the opposing brushes interdigitate 
because the distance between the opposing end-grafted 
beads is D = 17.5 < 2ho (ho denoting the unperturbed 
height of a single brush). 

Following previous works; 15 ' 21 ' 22 we quantify the 
amount of interdigitation via the overlap integral, I ov , 
defined as 



In 



p\(z) x p u (z)dz , 



(13) 



where p u and p\ are respectively the number densities of 
upper and lower brush layers, v mono =0.52 is the volume 
of a monomer, and A the surface area covered by the 
grafted beads. I ov follows from integrating the curves 
depicted in the insets of Fig. [T] Figure [7J shows this 
quantity as a function of the shear rate, v w /D, for good 
solvent conditions. The interdigitation is much higher for 
the larger grafting density, p g = 4.9/9* [Fig. [7(a)], as ex- 
pected from the more important stretching of the chains, 
and it becomes smaller with increasing shear rate and 
the progressive tilting of the chains. For the Langevin 
thermostat with 7 = 0.5 in all directions, this effect is 
more pronounced because the biasing of the flow profile 
increases the tilting of the brushes. A slight systematic 
difference can also be observed for the DPD thermostat 
with 7 = 0.5 where the overlap is systematically higher 
than in the other DPD cases. For 7 = 0.5 the tempera- 
ture is not properly conserved with the standard weight 
functions under good solvent conditions and the brush is 
additionally stretched. 

On the other hand, for p g = 1.2p*, the overlap is fairly 
constant for DPD and the Langevin thermostat with 
7^ = [Fig. [7(b)] over the whole interval of shear rates 
while for the standard Langevin thermostat (7^ ^ 0) the 
strong monotonous decrease is again related to the bias 
in shear direction. 

Figure [5(a) shows the shear stress a a = ((Fx) the 
mean force acting on the end-grafted beads in shear di- 
rection) times the inverse shear rate as a function of the 
constant relative wall velocity, v w , for good solvent con- 
ditions. As <j s D/v w reflects the "effective viscosity" of 
the polymeric system, the decrease of this quantity in- 
dicates a non-linear behavior, known as shear-thinning. 
The overall behavior of all cases is roughly the same. The 
Langevin thermostat approximately reproduces the DPD 
result after subtracting the constant shear stress given 
by the Langevin damping [Eq. (|12[) ] from the measured 
value. The Langevin thermostat with 7^ = quantita- 
tively agrees with the DPD case. The DPD thermostat 
was used with the standard weight functions for two dif- 
ferent values of the friction constant: 7 = 2 and 7 = 0.5. 
The latter value was also used in combination with con- 
stant weight functions. Only for 7 = 0.5, the standard 
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Figure 7: (Color online) Overlap integral versus shear rate 
for the two studied grafting densities. For DPD, we used the 
standard weight functions [Eq. (|10|l ] and the constant ones 
(uj = const). 



weight functions lead to some systematic differences for 
the shear stress at large wall velocities. We found out 
that this is due to the fact that under these conditions 
DPD fails to maintain the temperature at the desired 
value (T = 1.68). 

Figure[5(b) shows the effective viscosity for the smaller 
grafting density. The physical situation is now different 
as compared to the previous case: the opposing brush 
layers have a very small degree of interdigitation which is 
now independent of the wall velocity [see also Fig. [7(b)]. 

The linear response is observed for DPD and Langevin 
thermostats with 7^ = 0. For the standard Langevin 
thermostat (7^ ^ 0), however, the "effective viscosity" 
decreases and drops drastically for the largest wall ve- 
locity. This can be explained via the behavior presented 
in section Ull Al the tilting of the brush reduces the in- 
terdigitation of the brush layers not only because of the 
interaction among the brushes but also due to the strong 
interaction with the ficticious solvent that the Langevin 
thermostat inevitably implies. This is, however, an un- 
physical artifact in the simulations studied here. As a 
consequence, the density profiles depend on the parame- 
ters of the thermostat, which again is an unphysical ef- 
fect. Moreover, we emphasize that, even when the overall 
behavior of DPD and Langevin thermostats with j x = 
is similar, the absolute value of the shear stress is differ- 
ent. 
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Figure 8: (Color online) Comparison of "effective viscosity" for 
two opposing polymer brush layers under good solvent con- 
ditions for DPD and Langevin thermostats, using Langevin 
damping with 7^ 7^ and 7^ = . Different friction constants 
and two weight functions for DPD are considered. High and 
low grafting densities are presented in the upper and lower 
panel, respectively. 



Figure 9: (Color online) Comparison of normal pressure for 
two opposing polymer brush layers under good solvent condi- 
tions. DPD with different values of friction constants, 7, and 
weight functions are considered and compared to the standard 
Langevin thermostat (7 equal in all spatial directions) and the 
case 73; = 0. The upper panel shows a grafting density with 
strong interdigitation between the brushes, and the lower one 
presents a case with small brush-brush interdigitation. 



Figure [9] shows the normal stress as a function of 
shear rate for the two studied grafting densities. Fig- 
ure EJa) considers the higher grafting density, for which 
the brushes are strongly intcrdigitatcd and slightly com- 
pressed (under good solvent conditions the mean force 
between the layers is repulsive). A decreasing normal 
stress as a function of shear velocity is found, except for 
the case of DPD with 7 = 0.5 and the standard weight 
functions. As mentioned above, for this case DPD does 
not properly conserve temperature and a slight heat up 
of the system is observed, which in turn produces a fur- 
ther increase in the steady state stretching of the brush 
with a concomitant increase in the normal repulsion of 
the brush layers. For all the other cases, the decrease of 
normal pressure upon increasing velocity is produced by 
the progressive tilting of the chains and the decrease of 
interdigitation, already observed in the behavior of the 
overlap integral [Fig. EJa)]. 

Figure [5{b) shows the normal stress for the smaller 
grafting density. Here, the brush is so dilute that the in- 
teraction between the brush layers is almost negligible. In 
this case, there is a mean attraction between the layers, 
due to the wall interaction with each bead [see Eq. Q]. 
For DPD (7 = 2.5) and the Langevin thermostat with 



7a: = 0, the structure of the brush is quite similar result- 
ing in a similar behavior of the normal force. The small 
interdigitation leads to a very small change of the incli- 
nation angle giving rise to a very weak dependence of the 
normal stress on the wall velocity. This behavior agrees 
with the approximately constant behavior of the overlap 
integral shown in Fig. [7f b) and the linear response ob- 
served in the effective viscosity. A different behavior is 
found for the Langevin case with 7 = 0.5 in all directions. 
Under these conditions, the unphysical enhancement of 
the chain inclination leads to a much larger interaction 
between wall and monomers, which is mainly attractive 
for typical bead positions. 

Normal and shear stresses were also investigated under 
poor solvent conditions. We used the Langevin thermo- 
stat with 7 = 0.5 perpendicular to the shear direction 
and 73; = 0. For the DPD thermostat with 7 = 0.5, the 
temperature is conserved under poor solvent conditions 
unlike in the good solvent case. As already shown in ta- 
ble HI the larger cut-off radius for poor solvent conditions 
improves the DPD efficiency and keeps T constant even 
for the standard weight functions. 

Figure HOf a) shows a very similar behavior of both 
thermostats concerning the effective viscosity, I ov , and 
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Figure 10: (Color online) Comparison of shear and normal 
stress for DPD and Langevin thermostat with -y x = (poor 
solvent conditions). Both thermostats give a similar result. 
The inset shows the overlap integral I ov as defined in Eq. (|13fl . 



the normal stress. Within the considered regime of shear 
rates, both thermostat implementations lead to equiva- 
lent results. As compared with the good solvent case, the 
normal stress at the wall is negative, i.e. there is a mean 
attraction between upper and lower brush layers. This is 
due to the mean attraction among beads for this model. 

If this equivalence between DPD and LGV with 7 T = 
holds also for stronger out-of-equilibrium situations, i.e at 
higher shear velocities, remains to be investigated thor- 
oughly. 



The role of weight functions in the DPD 
thermostat 



In this section we consider, in more detail, the ability 
of the DPD thermostat to conserve temperature in non- 
equilibrium simulations for different weight functions. 
The system under study is a polymeric liquid, formed 
by 10-bead chains, confined by two polymer brushes of 
identical chains.— The brush and melt density profiles 
across the perpendicular directions correspond to that of 
Fig. [TJc). We imposed either a Poiseuille flow by means 
of a constant external volume force or a linear Couette 
flow by moving the walls at constant relative velocity. 

In Fig. 111! the violation of temperature conservation 
is shown as a function of external force, f x , for Poiseuille 
flow using the standard DPD weight functions. The tern- 
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Figure 11: (Color online) Poiseuille flow for different exter- 
nal forces, f x . Upper panel: Temperature, as obtained from 
the mean square velocity perpendicular to the shear direc- 
tion. Lower panel: Velocity profile across the polymeric liq- 
uid. Only for the two smallest external forces DPD with the 
standard weight functions [Eq. (|10[) ] maintains the temper- 
ature at the desired value. The simulations were performed 
under poor solvent conditions. 



perature, liquid number density, and brush grafting den- 
sity were respectively set to T = 1.68, p m = 0.61 and 
p s = 5.5p* (= 0.77a- 2 ), with p* = l/nR 2 g and R CJ = 1.50. 
The typical Poiseuille velocity profile across the gap is 
shown in Fig. lllf b). while the temperature profile - mea- 
sured by the mean square velocity in the direction per- 
pendicular to the flow - is presented in Fig. fTTI a). Tem- 
perature is conserved only for the two smallest forces. 
In the remaining cases, the temperature increases in the 
region of large velocity gradients. These cases show ex- 
amples in which the DPD thermostat fails to maintain 
the desired temperature even under poor solvent condi- 
tions. 

We also studied the temperature conservation for Cou- 
ette flows using the weight functions considered in table[l] 
In Fig. [r2T a) the temperature profiles for a shear velocity 
of v w = 3 arc shown. 

The standard weight functions clearly fail to keep tem- 
perature constant and lead to quadratic temperature pro- 
files with a maximum in the middle of the gap. In con- 
trast to the Poiseuille flow this is not related to the 
velocity gradient, which is constant across the gap for 
Couette flows. We attribute the resulting temperature 
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Figure 12: (Color online) Comparison of different DPD weight 
functions for a Couette flow (poor solvent conditions) . In the 
upper panel, v w — 3 is considered: while the standard choice 
of DPD weight functions is not able to conserve temperature, 
a choice of lor = (1 — 7~) 1//2 improves temperature conserva- 
tion, and lor = Q(r c — r) conserves temperature at the chosen 
value (indicated by the dashed line). At v w = 8 (lower panel 
), temperature is not conserved at the desired value, and the 
behavior across the film is very different for the considered 
weight functions (see text). 



profiles to the density distribution of monomers, which 
is enhanced close to the brush coated walls giving rise 
to a local improvement of the efficiency of the ther- 
mostat in this region. Another choice of weight func- 
tions [u; R = VuP = (1 — r/rc) 1 / 2 ] gives a better re- 
sult although it also fails to conserve temperature. The 
constant- weight functions (lo r = V uP = 1) arc more 
efficient and conserve temperature. 

Figure [LWb) shows the temperature profile for a shear 
velocity of v w = 8. Under this condition, the thermostat 
is not able to conserve temperature because for reason- 
able values of the friction constant energy cannot be dis- 
sipated as fast as it is plugged into the system. The solid 
line, corresponding to the standard DPD weight func- 
tions, shows a similar behavior as in the previous case: 
temperature is conserved more efficiently in the regions 
of higher density. For constant weight functions, temper- 
ature is fairly constant all across the film but differs from 
the desired value (indicated by a dashed line). 

>From these examples of strong out-of-equilibrium 
simulations, we conclude that the choice of weight func- 
tions can make an important difference. For example, 



the standard weight functions might not be able to con- 
serve temperature even for a physically meaningful case. 
Constant weight functions seem to be a good alterna- 
tive and are even more efficient computationally, as al- 
ready noted in a previous study^ As a general guideline, 
it should be checked that the physical conditions fulfill 
the relation 7 = A^tp7dpd/w, where 7 is the shear stress 
imposed to the system. In any case, for short range po- 
tentials (as Lennard- Jones and good solvent conditions) 
or strong out-of-equilibrium simulations, the temperature 
profile across the sample, as shown in Figure [T^] can give 
insight on the ability of the thermostat for keeping the 
temperature constant. 



IV. CONCLUSIONS 

In this work we tested and compared commonly used 
implementations of Langevin and DPD thermostats for 
different polymeric systems. Equilibrium, transient and 
steady state conditions were considered for the study of 
various reference systems, such as polymer brush bear- 
ing surfaces or brush-melt interfaces. We utilized a well 
studied coarse-grained bead-spring polymer model. By 
varying the cut-off in the interaction potential we mim- 
icked good and poor solvent conditions. We quantified 
the relative strength of the thermostats in a wide range 
of friction constants, 7, and found that the strength of 
the Langevin thermostat is much larger than DPD with 
standard weight functions for similar values of 7. The 
simulation of the transient state of a polymer brush layer 
driven to constant velocity from rest illustrates the known 
weaknesses of the Langevin thermostat - lack of mo- 
mentum conservation, screening of hydrodynamic inter- 
actions, and violation of Galilean invariance - and how 
these are avoided by the DPD thermostat which con- 
serves local momentum. When applied in shear direction, 
the Langevin thermostat biases the velocity profile. The 
common workaround of switching-off the Langevin ther- 
mostat in the non-equilibrium direction was analyzed for 
different systems. We found that, in most cases, the lat- 
ter behaves similar to the DPD scheme but care must be 
exerted when the system is strongly driven out of equi- 
librium. 

We furthermore quantified the differences between var- 
ious forms of weight functions of the DPD thermostat, 
which can be chosen freely, provided that the weights 
for random and dissipativc forces obey relation It 
is important to note that most of the previous works us- 
ing DPD utilized a standard form [Eq. (|10p]. which was 
originally intended to be used in conjunction with "soft" 
potentials. When the DPD thermostat is applied to typ- 
ical coarse-grained potentials, the "hard" nature of the 
conservative potentials prevents one from using a very 
large time step, and therefore not much is gained from 
"smooth" thermostat forces. Moreover, we found that the 
typical weights can be regarded as adequate for equilib- 
rium and slightly out-of-equilibrium conditions but they 
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fail to conserve temperature for medium and strong driv- 
ing forces. We tested this for both, Poiseuille and Couette 
flows, of a polymer melt confined between two polymer 
brush layers. We found quantitatively that taking con- 
stant weight functions is both computationally faster and 
yields a thermostat that is more suitable for strong out- 
of-equilibrium situations in which a large amount of heat 
per unit time is produced. 

It is important to note that none of the methods dis- 
cussed in the present paper is suitable to fully account for 
hydrodynamic effects at arbitrary densities. For instance, 
the coupling between the monomer density distribution 
and the external flow profile is not described. To achieve 
this, other methods, e.g., the self-consistent solution of 
the Brinkman equation^! have to be applied. In a future 
study we plan to investigate the systems considered here 
using explicit solvent molecules. This should put us into 
the position to understand the importance of the effects 
delineated above. 

Finally, our results clearly show that great care is 
needed in non-equilibrium MD simulations of soft matter 



systems in order to ensure that the simulations are free 
of artifacts due to an inappropriate choice of the thermo- 
stat. 
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